2022-03-13: Myelin phagocytosis

Data

Macrophages derived from Control (CT) and PD individuals were exposed to myelin prepared from CT individuals.

df = read_xlsx("ctrl_myelin-macrophages.xlsx")

df

df %<>%
  pivot_longer(-TIME, values_to = "INT", names_to = "ID") %>%
  separate(ID, sep = 2, into = "STATUS", remove = FALSE) %>%
  select(ID, STATUS, TIME, INT) %>%
  arrange(ID, TIME)

df

ggplot(data = df, mapping = aes(x = TIME, y = INT, color = STATUS)) + geom_smooth(method = "loess", formula="y ~ x")

ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = STATUS)) +
  geom_point() +
  geom_line()

Linear regression forcing intercept to zero

ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = STATUS)) + geom_smooth(method = "lm", formula="y ~ 0 + x", se = FALSE)

m1.fit = lmer(INT ~ 0 + TIME*STATUS + (0 + TIME | ID), data=df)

Warning in as_lmerModLT(model, devfun) :
  Model may not have converged with 1 eigenvalue close to zero: 1.3e-11

summary(m1.fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ 0 + TIME * STATUS + (0 + TIME | ID)
   Data: df

REML criterion at convergence: 14800

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.152 -0.443 -0.061  0.439  7.594 

Random effects:
 Groups   Name Variance        Std.Dev.
 ID       TIME  13955019515874  3735642
 Residual      120010088567134 10954912
Number of obs: 420, groups:  ID, 20

Fixed effects:
                                Estimate                 Std. Error                         df t value             Pr(>|t|)    
TIME                           6081194.1                  1132610.7                       18.3    5.37       0.000039822857 ***
STATUSCT                      17763784.6                  1391548.9 114154268801901680132096.0   12.77 < 0.0000000000000002 ***
STATUSPD                      10266197.3                  1538415.2  79120202192516307484672.0    6.67       0.000000000025 ***
TIME:STATUSPD                 -2356761.2                  1688396.3                       18.3   -1.40                 0.18    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            TIME   STATUSC STATUSP
STATUSCT    -0.090                
STATUSPD     0.000  0.000         
TIME:STATUS -0.671  0.060  -0.067 

m1.res = simulateResiduals(fittedModel = m1.fit)
plot(m1.res)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

Linear regression forcing intercept to zero and up to 10 hours

m2.fit = lmer(INT ~ 0 + TIME*STATUS + (0 + TIME | ID), data=df %>% filter(TIME <= 10))

Warning in as_lmerModLT(model, devfun) :
  Model may not have converged with 1 eigenvalue close to zero: 1.2e-11

summary(m2.fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ 0 + TIME * STATUS + (0 + TIME | ID)
   Data: df %>% filter(TIME <= 10)

REML criterion at convergence: 7589

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.268 -0.413 -0.039  0.234 10.734 

Random effects:
 Groups   Name Variance       Std.Dev.
 ID       TIME 18017163807632 4244663 
 Residual      64442444229858 8027605 
Number of obs: 220, groups:  ID, 20

Fixed effects:
                                Estimate                 Std. Error                         df t value  Pr(>|t|)    
TIME                           8716594.4                  1300454.6                       18.8    6.70 0.0000022 ***
STATUSCT                       5488175.9                  1365297.6   1883035153459291881472.0    4.02 0.0000583 ***
STATUSPD                       4298811.9                  1509393.3 420829932100792145149952.0    2.85    0.0044 ** 
TIME:STATUSPD                 -3747708.4                  1938603.3                       18.8   -1.93    0.0684 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            TIME   STATUSC STATUSP
STATUSCT    -0.150                
STATUSPD     0.000  0.000         
TIME:STATUS -0.671  0.101  -0.111 

m2.res = simulateResiduals(fittedModel = m2.fit)
plot(m2.res)

Linear regression not forcing intercept to zero

ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = STATUS)) + geom_smooth(method = "lm", formula="y ~ x", se = FALSE)

m3.fit = lmer(INT ~ TIME*STATUS + (TIME | ID), data=df)

Warning in as_lmerModLT(model, devfun) :
  Model may not have converged with 1 eigenvalue close to zero: 1.4e-11

summary(m3.fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ TIME * STATUS + (TIME | ID)
   Data: df

REML criterion at convergence: 14777

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.257 -0.338 -0.020  0.384  7.875 

Random effects:
 Groups   Name        Variance        Std.Dev. Corr
 ID       (Intercept)  46033779618314  6784820     
          TIME         12199710583965  3492808 0.44
 Residual             108281639341729 10405846     
Number of obs: 420, groups:  ID, 20

Fixed effects:
              Estimate Std. Error       df t value   Pr(>|t|)    
(Intercept)   17763785    2435581       18    7.29 0.00000089 ***
TIME           6081194    1059174       18    5.74 0.00001922 ***
STATUSPD      -7497587    3630750       18   -2.07      0.054 .  
TIME:STATUSPD -2356761    1578923       18   -1.49      0.153    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) TIME   STATUS
TIME         0.315              
STATUSPD    -0.671 -0.211       
TIME:STATUS -0.211 -0.671  0.315

m3.res = simulateResiduals(fittedModel = m3.fit)
plot(m3.res)

Linear regression not forcing intercept to zero and up to 10 hours

m4.fit = lmer(INT ~ TIME*STATUS + (TIME | ID), data=df %>% filter(TIME <= 10))

Warning in as_lmerModLT(model, devfun) :
  Model may not have converged with 1 eigenvalue close to zero: 1.2e-11

summary(m4.fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ TIME * STATUS + (TIME | ID)
   Data: df %>% filter(TIME <= 10)

REML criterion at convergence: 7587

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.642 -0.387 -0.020  0.213 10.810 

Random effects:
 Groups   Name        Variance       Std.Dev. Corr
 ID       (Intercept)  9830589540884 3135377      
          TIME        17362487055394 4166832  0.12
 Residual             61633698670715 7850713      
Number of obs: 220, groups:  ID, 20

Fixed effects:
              Estimate Std. Error       df t value  Pr(>|t|)    
(Intercept)    5488176    1635996       18    3.35    0.0035 ** 
TIME           8716594    1276458       18    6.83 0.0000022 ***
STATUSPD      -1189364    2438799       18   -0.49    0.6317    
TIME:STATUSPD -3747708    1902831       18   -1.97    0.0645 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) TIME   STATUS
TIME        -0.052              
STATUSPD    -0.671  0.035       
TIME:STATUS  0.035 -0.671 -0.052

m4.res = simulateResiduals(fittedModel = m4.fit)
plot(m4.res)

Two-factor repeated-measures ANOVA

ggplot(data = df, mapping = aes(x = factor(TIME), y = INT, color = STATUS)) +
  geom_boxplot()

summary(aov(INT ~ STATUS*factor(TIME) + Error(ID), data = df))

Error: ID
          Df             Sum Sq            Mean Sq F value Pr(>F)  
STATUS     1 100316595603111792 100316595603111792    3.23  0.089 .
Residuals 18 558649433744305344  31036079652461408                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
                     Df             Sum Sq           Mean Sq F value              Pr(>F)    
factor(TIME)         20 405726546818305024 20286327340915252   37.89 <0.0000000000000002 ***
STATUS:factor(TIME)  20  23064109848552772  1153205492427638    2.15               0.003 ** 
Residuals           360 192757182833994944   535436618983319                                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(lmer(INT ~ STATUS*factor(TIME) + (1 | ID), data=df))

Warning in as_lmerModLT(model, devfun) :
  Model may not have converged with 1 eigenvalue close to zero: 2.7e-12

Type III Analysis of Variance Table with Satterthwaite's method
                                Sum Sq           Mean Sq NumDF DenDF F value              Pr(>F)    
STATUS                1729285478721196  1729285478721196     1    18    3.23               0.089 .  
factor(TIME)        382825994343270720 19141299717163536    20   Inf   35.75 <0.0000000000000002 ***
STATUS:factor(TIME)  23064109848552360  1153205492427618    20   Inf    2.15               0.002 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(lmer(INT ~ STATUS*factor(TIME) + (1 | ID), data=df))

Warning in as_lmerModLT(model, devfun) :
  Model may not have converged with 1 eigenvalue close to zero: 2.7e-12

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ STATUS * factor(TIME) + (1 | ID)
   Data: df

REML criterion at convergence: 14062

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.528 -0.571 -0.004  0.487  3.333 

Random effects:
 Groups   Name        Variance         Std.Dev.
 ID       (Intercept) 1452411626813050 38110519
 Residual              535436618009765 23139503
Number of obs: 420, groups:  ID, 20

Fixed effects:
                                       Estimate              Std. Error                      df t value             Pr(>|t|)    
(Intercept)                             31620.1              13442971.3                    32.6    0.00              0.99814    
STATUSPD                               -21728.0              20039598.4                    32.6    0.00              0.99914    
factor(TIME)1                         8781401.4               9866717.4  42304711666901590016.0    0.89              0.37346    
factor(TIME)2                        23768179.7               9866717.4  41019441770847150080.0    2.41              0.01600 *  
factor(TIME)3                        34487435.5               9866717.4  43569307823568625664.0    3.50              0.00047 ***
factor(TIME)4                        50513832.0               9866717.4  38613191227455332352.0    5.12     0.00000030615385 ***
factor(TIME)5                        51811985.3               9866717.4  45441836676151836672.0    5.25     0.00000015112156 ***
factor(TIME)6                        58989519.0               9866717.4  40347422911320178688.0    5.98     0.00000000225013 ***
factor(TIME)7                        66618580.2               9866717.4  40132082484553932800.0    6.75     0.00000000001460 ***
factor(TIME)8                        73883423.8               9866717.4  38849758612996620288.0    7.49     0.00000000000007 ***
factor(TIME)9                        82034966.7               9866717.4  39666229042097963008.0    8.31 < 0.0000000000000002 ***
factor(TIME)10                       88545483.7               9866717.4  38580371851014479872.0    8.97 < 0.0000000000000002 ***
factor(TIME)11                       92236016.1               9866717.4  43039447898024665088.0    9.35 < 0.0000000000000002 ***
factor(TIME)12                       95729228.8               9866717.4  38863621046500278272.0    9.70 < 0.0000000000000002 ***
factor(TIME)13                      102929904.9               9866717.4  40508903641680633856.0   10.43 < 0.0000000000000002 ***
factor(TIME)14                      106728537.7               9866717.4  41475839895920320512.0   10.82 < 0.0000000000000002 ***
factor(TIME)15                      111739025.5               9866717.4  43975724117093941248.0   11.32 < 0.0000000000000002 ***
factor(TIME)16                      112668644.3               9866717.4  41299628101740830720.0   11.42 < 0.0000000000000002 ***
factor(TIME)17                      117035519.1               9866717.4  42709417196245622784.0   11.86 < 0.0000000000000002 ***
factor(TIME)18                      119225475.0               9866717.4  41683750762701217792.0   12.08 < 0.0000000000000002 ***
factor(TIME)19                      123280234.6               9866717.4  40807518497940914176.0   12.49 < 0.0000000000000002 ***
factor(TIME)20                      128418824.0               9866717.4  42496879823159566336.0   13.02 < 0.0000000000000002 ***
STATUSPD:factor(TIME)1               -2192325.3              14708433.9 302803913031962787840.0   -0.15              0.88151    
STATUSPD:factor(TIME)2               -7940390.1              14708433.9 263360823376362176512.0   -0.54              0.58930    
STATUSPD:factor(TIME)3              -11572207.9              14708433.9 315759423867243003904.0   -0.79              0.43141    
STATUSPD:factor(TIME)4              -23749022.1              14708433.9 286506603967034687488.0   -1.61              0.10639    
STATUSPD:factor(TIME)5              -20600772.6              14708433.9 311933493985617313792.0   -1.40              0.16133    
STATUSPD:factor(TIME)6              -23914955.2              14708433.9 321361150475973558272.0   -1.63              0.10396    
STATUSPD:factor(TIME)7              -27324733.9              14708433.9 243030239043486711808.0   -1.86              0.06320 .  
STATUSPD:factor(TIME)8              -30697005.5              14708433.9 242991815025575854080.0   -2.09              0.03689 *  
STATUSPD:factor(TIME)9              -33806339.4              14708433.9 279563481426973523968.0   -2.30              0.02154 *  
STATUSPD:factor(TIME)10             -37170207.8              14708433.9 241790702101374074880.0   -2.53              0.01150 *  
STATUSPD:factor(TIME)11             -39111121.5              14708433.9 289315740243068125184.0   -2.66              0.00784 ** 
STATUSPD:factor(TIME)12             -37386024.8              14708433.9 235730788509275062272.0   -2.54              0.01103 *  
STATUSPD:factor(TIME)13             -40031013.0              14708433.9 242992509354682646528.0   -2.72              0.00650 ** 
STATUSPD:factor(TIME)14             -38876187.6              14708433.9 270572760748935938048.0   -2.64              0.00821 ** 
STATUSPD:factor(TIME)15             -42234173.8              14708433.9 313696440606546329600.0   -2.87              0.00409 ** 
STATUSPD:factor(TIME)16             -43377542.1              14708433.9 273472501187227844608.0   -2.95              0.00319 ** 
STATUSPD:factor(TIME)17             -46923491.4              14708433.9 306408864699436498944.0   -3.19              0.00142 ** 
STATUSPD:factor(TIME)18             -47063593.3              14708433.9 292368502840325636096.0   -3.20              0.00138 ** 
STATUSPD:factor(TIME)19             -48010183.6              14708433.9 271948042683244412928.0   -3.26              0.00110 ** 
STATUSPD:factor(TIME)20             -49931609.0              14708433.9 284229729013514764288.0   -3.39              0.00069 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 42 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

Model curve using aomisc package

m = df %>%
  nest(data = -c(ID, STATUS)) %>%
  mutate(
    model = map(.x = data, ~ drm(.x$INT ~ .x$TIME, fct = DRC.negExp())),
    tidied = map(model, tidy)
  )

m

tt = m %>% 
  unnest(tidied) %>%
  pivot_wider(id_cols = c(ID, STATUS), names_from = term, values_from = estimate) %>%
  rename(PLATEAU = a, RATE = c)

tt

t.test(PLATEAU ~ STATUS, data = tt)

    Welch Two Sample t-test

data:  PLATEAU by STATUS
t = 0.39, df = 15, p-value = 0.7
alternative hypothesis: true difference in means between group CT and group PD is not equal to 0
95 percent confidence interval:
 -130349196  188986685
sample estimates:
mean in group CT mean in group PD 
       175235687        145916942 

ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = ID)) +
  geom_point() +
  geom_line()

---

See also:

• https://stats.oarc.ucla.edu/r/seminars/repeated-measures-analysis-with-r/
• https://doi.org/10.1213%2FANE.0000000000003511
• https://cran.r-project.org/web/packages/growthrates/vignettes/Introduction.html
• http://dpmartin42.github.io/posts/Piecewise-growth
• https://doi.org/10.1016/j.buildenv.2021.108349
• https://doi.org/10.3758/s13423-017-1281-0